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Abstract: Recently, 3D image fusion reconstruction using a FDK algorithm along three- 
orthogonal circular isocentric orbits has been proposed. On the other hand, we know that 3D 
image reconstruction based on three-orthogonal circular isocentric orbits is sufficient in the 
sense of Tuy data sufficiency condition. Therefore the datum obtained from three- orthogonal 
circular isocentric orbits can derive an exact reconstruction algorithm. In this paper, an exact 
weighted-FBP algorithm with three-orthogonal circular isocentric orbits is derived by means 
of Katsevich's equations of filtering lines based on a circle trajectory and a modified weighted 
form of Tuy's reconstruction scheme. 
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1. Introduction 

The cone-beam scanning configuration with a circular trajectory remains one of the most popular 
scanning configuration and has been widely employed for data acquisition in 3-D X-ray computed to- 
mography (CT), because it allows for operating at a high rotating speed due to its symmetry, avoiding 
the need to axially translate the patient, such as in helical or step-and-shoot CT [1, 2]. Unfortunately, a 
circular trajectory does not satisfy Tuy's data sufficiency condition [3] which requires every plane that 
passes through the reconstruction region (ROI) must also cut the trajectory at least once and therefore 
only approximate reconstructions are possible. In the pursuit of the data sufficiency condition, various 
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scanning trajectories have been proposed, such as circle and line [4, 5], circle plus arc [6, 7], double 
orthogonal circles [8, 9] and dual ellipses [10]. 

Recently, the FDK algorithm [11] has been implemented for a cone-beam vertex trajectory consisting 
of three-orthogonal isocentric circles [12]. On the other hand, we know that 3D image reconstruction 
based on three-orthogonal circular isocentric orbits is sufficient in the sense of Tuy data sufficiency 
condition. Therefore the datum obtained from three-orthogonal circular isocentric orbits can derive 
an exact reconstruction algorithm. In this paper, another theoretically exact and general inversion for- 
mula which was proposed in Katsevich [13] will be implemented for three-orthogonal circular isocentric 
orbits. Compared with the aforementioned algorithms, the distinctive features of Katsevich's algorithm 
can be summarized as the choice of a more general weight and a novel way to dealing with discontinuous 
in weight. Furthermore, Katsevich's algorithm has been implemented for various scanning trajectories, 
such as circle and line [14, 15], circle plus arc [16], several circular segment [17], two orthogonal circles 
[13] and helix [18, 19]. 

In contrast to a singular circular trajectory, there are several advantages in using weighted reconstruc- 
tion with the three-orthogonal-circular trajectory. First, for the reconstruction noise due to the algorithm 
is 'trajectory' dependent. Using three-orthogonal geometry with weighting function, the noise is sup- 
pressed. Second, the three-orthogonal setting yields better image quality in comparison with classical 
scheme for larger cone-angles. Last, the scan method for three-orthogonal geometry can be implemented 
with minimal requirements and can provide a better 3D reconstruction for small animals or objects. At 
the same time, compared to the two-circular trajectory, artifacts of the boundaries of the ROI where 
artifacts are likely to appear are suppressed by the additional filtering lines and according weighting 
functions. 

In our implementation, a new weighting scheme is developed so that all measurements are used 
by accurately averaging over multiply measured projections of the three-orthogonal-circular scanning. 
Moreover, this new weighted reconstruction is differs from the weighted FDK [11] reconstruction with 
the three-orthogonal-circular trajectory [12] considerably, in that it is a theoretically exact formulation of 
a general shift invariant filtered back-projection (FBP) reconstruction framework. In this paper, we show 
the derivation of our cone-beam FBP reconstruction algorithm starting with Tuy's inversion formula. In 
other words, we discuss how to derive a FBP cone-beam reconstruction formula from the Tuy's classical 
inversion formula. To construct the reconstruction formula, several operators are needed. First, using 
some properties of the cone-beam transform, the classical Tuy's inversion formula can be written as (the 
first-order derivative of the Radon transform in IR 3 ) Grangeat's inversion formula [20] in a weighted 
form. Then, using the methodology in[13], we can get the final inversion formula. The features of our 
inversion formula can be summarized as follow: using some properties of the cone-beam transform to 
derive the inversion formula, a proper choice of the weighting function, deriving equations of the filtering 
lines and describing geometric properties of filtering lines in the planar detector plane. 

The organization of this paper is as follows. In Section 2, we derive Katsevich's inversion formulation 
by Tuy's weighted form. In Section 3, we derive equations of the filtering planes and filtering lines. In 
Section 4, we show the geometric properties of filtering lines in the planar detector plane and according 
weighting functions. 
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2. Weighting Function for Tuy's Formula 

Let a trajectory C be a differentiable curve in IR 3 described by y(s), s E K. The object density function 
is /(x), where a; is a vector in the (yi, y 2 , Z/3) coordinate system, and /(x) is an infinitely differentiable 
real integrable function with a compact support f2 C R 3 \C. The modified cone-beam projection of f(x) 
along the direction of a/||a|| at the focal point location y(s) is defined as: 

Df(y(s), a) = f(y + ta)dt =-±-Df(y(s), ^-), (y(s), a) E C x S 2 . (1) 

7o IMI IMI 

Then, let £>/ be extended from C x § 2 to I 3 x R 3 as: 

D y f(z) = Df(y,z),y,zeR 3 . (2) 
First consider the three-dimensional Fourier transform of D y f(z) for a fixed /(x) as: 

DJ(a)= [ Df(y,z)e- 2mz -°dz. (3) 
Jr 3 

It can be easily verified that: 

D y ~f(ta) = ^D y ~f( ( T),for t > 0. (4) 

In the following, we will show how to derive a new reconstruction formula from Tuy's reconstruction 
scheme in a weighted form. 

Theorem 1. [Tuy's weighted form] Let C be a curve in R 3 parameterized by a piecewise differentiable 
function y(s), s el. For a fixed re G f2, we suppose that there exists a weighting function n x (s,a) : 
R x § 2 — > R such that ra x (s, a) is integrable with respect to the second variable a E E> 2 for each sei, 
and the set: 

R(x, a) = {s E R \x ■ a = y{s) ■ a, y'(s) ■ o ± 0} (5) 
is non-empty and finite for almost all a E § 2 . Then: 

provided that n x (s, a) fulfilled the completeness condition: 

n x (s,a) = l,a.e.ina E § 2 . (7) 

sS-R(x,ct) 

This inversion formula also tells us Tuy's data sufficiency conditions for an accurate reconstruction 
of a ROI from cone -beam projections. These conditions are: a ■ y(s) = a ■ x and a ■ y'(s) 7^ 0. In the 
following, we will show how to derive a FBP reconstruction formula from Tuy's formula and show how 
Tuy's data sufficiency and nonended condition can be relaxed. In a first step we try to use Equation (1) 
to simplify Equation (4) as follows: 
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iD y f(a) = i [ Df(y(q), z)e~ 2 ^dz = i I [°° -Df(y(q), 9)e~ 2 ^r 2 drd9, (8) 
Jr3 Js 2 Jo r 

where r is the spherical radical coordinate, 9 is a three dimensional unit vector such that z = r6 and 
|| #|| = 1. Next we introduce the quantity: 

iDJ(a) = lf Df(y(q),6)d9 T ire~ 2 ™ e -°dr + l [ Df(y(q),9)d9 H ' i\r\e~ 2 ^dr. (9) 

1 JS 2 J~oo 1 JS 2 J-oo 

In Equation (7) the first term is real and add, and the second term is imaginary and even if we take 
into account the fact that the factor y'(s) ■ a in the inversion formula is also odd in a, we conclude that 
only the first term contributes in Tuy's inversion formula. The second term will vanish when we perform 
the integration over a. Thus we have: 

iDJ(a) = ^ / Df(y(q),8)6'(6 ■ a)d9. (10) 
Thus, we rederived the weighted form of Tuy's inversion formula as follows: 

Therefore, there is a derivative of delta function in the integration over 9. We can perform integration 
by part for the integration over 9 in Equation (10). We thus obtain the directional derivative along the 
direction a of cone-beam data Df(y(q), 9) with respect to unit vector 9. Therefore Equation (10) can be 
written as follows: 

= E ¥y^I"(/ Ve,*Df(y( q ),6)8(e-(j))\ q=s ded<T. (12) 

Using Grangeat's formula and change of variables p — > q defined by p = y(q) ■ a, we obtain: 



Ve,aDf{y{q),9)de\ \ q=s = a)\ p=y{s) . a=x . a , (13) 



i d r r _ „ „ i d 2 

y'(s) ■ a dq 

where f(p, a) is the 3D Radon transform of /. Modifications from Tuy's original inversion (6) to Equa- 
tion^) with change of variables have been important progress. Accounting for the numerical imple- 
mentation of formula Equation (11) is still complicated. We can change the two integrals over the unit 
vector 9 = (— sin 7 cosy?, sin 7 cosy?, cos 7) and a = (cosy?, sin y?, 0) into three single integrals over 
parameters s, ip and 7. Using the methodology in [13], we get: 



f{x) = J] C ^ {sgn{y'{s) ■ a)n x {s, a)}dip 

x J 0 2W ^ jD /(^(^)' cos ^( s ' x ) + sm 7e*(s,y))| (?=s ^ 
where (3(s,x) = {x- y{s))/\\x - y(s)\\ = (0,0, l),e x (s,<p) = (3(s,x) x a x (s,<p), denote 

k x (s, ip) = sgn(y'(s) ■ a)n x (s, ip) 

w m (s, x) = ]ha(k x (s, ip m + e) - k x (s, <p m - e)). 
e^o 



(14) 



(15) 
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The values of w m (s, x) are determined by the definition of the weighting function and signum function 
near the discontinuous points (f m 's. Furthermore, in the FBP reconstruction scheme, the unit vector 
&x(s, tp m ) is the normal vector of filtering plane and the vector e x (s, (p m ) denotes the direction of filtering 
lines which pass through the object point x and the source point y(s). After substituting Equation (14) 
into formula Equation (13), we obtain: 

1 f w m(s,x) f 2jT d , dr 

f( x ) : =-^ / y.] r\\ x / ^Df(y(q),cosrP(s,x) + smre x (s,ip m ))\ q =8- — as. (16) 

°^ J I \ x ~ U\ s )\ io smr 

The final inversion formula is independent of the specific geometrical shape of the image object. 
Rather, it is determined by the scanning geometry. Thus, we need study the specific features of filtering 
lines and weighting functions. In the following, the inversion formula will be implemented for a cone- 
beam vertex trajectory consisting of three-orthogonal to each other circular. Furthermore, the filtering 
lines and the weighting functions will be given in detail. 

3. Equation of Filtering Lines for the Three-Orthogonal Circular Scanning 

We propose the use of three independent, orthogonal to each other, circular isocentic scan-paths for 
X-ray projection acquisitions. Let a trajectory C = G\ U C 2 U C 3 , as shown in Figure 1. 

Figure 1. Geometries for three-orthogonal circular trajectory and the planar-detector plane. 




where: 

Ci = {y(s) G R 3 : yi (s) = Rcoss,y 2 (s) = Rsms,y 3 (s) = 0, s G [0,2tt]}; (17) 

C 2 = {y(s) G -R 3 : y 1 (s) = Rcoss,y 2 (s) = 0,y 3 (s) = Rsms,s£ (2n,An]}; (18) 

C 3 = {y(s) G R 3 : 2/i (s) = 0,y 2 (s) = Rcoss,y 3 (s) =Rsins,s e (47r,67r]}. (19) 

Let f(x) is the density function to be constructed. Assume that the function is smooth and vanishes 
outside the ball: 



tt = {x = (y a , y 2 , y 3 )\y{ + y\ + y\ < r}, 0 < r < R, 



(20) 
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where r is the radius of the subject ball and R is the radius of the scanning circular. We assume that the 
physical detector is a planar-detector, which is denoted as DP(s), where s is the parameter of source 
point. Furthermore, let the planar-detector is at distance 2R from the source, as shown in Figure 1. 
Denote by (u, v) the horizontal and vertical coordinates of the planar-detector plane, where u is parallel 
to the tangent vector of source, v is the normal vector of scanning trajectory and the origin is at the 
projection of the source point y(s) onto the detector plane. 

Therefore, in the (yi, yi-, yz) coordinate system, any point on the detector plane can be characterized 
completely by use of (u, v). It can be shown that: 

yx — u sin s — R cos s, y<i = — w cos s — R sin s, y^ = v, s G [0, 27r]; (21) 
y\ = m sin s — R cos s,y2 — v, y% = —u cos s — R sin s, s G (2ir, 47r] ; (22) 
yi = v, y-i — u sin s — R cos s,y% = —u cos s — R sin s, s G (47T, 6n] . (23) 

Figure 2. The filtering plane is shown, which is tangent to one of circular trajectories and is 
perpendicular to a x . 




Figure 2 shows a filtering plane R(x, <r)through the y(s) and is tangent to other trajectory at point 
y(X). Without loss of generality we consider the source y(s) G C\ (the case y(s) G 6*2,3 is treated 
in a similar fashion). The point of tangency is given by either y(X) = (i?cos A, 0, i?sin A) or y(X) = 
(0, R cos A, R sin A) . In the first case, the normal vector to the tangent plane is given by: 

a x = (sin s cos A, 1 — cos s cos A, sin s sin A) . (24) 
The equation of tangent plane, which is a filter plane, is given by: 



R(x, a) = {x \(x — y(s)) ■ a = 0} 



(25) 
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Filtering lines on the detector are obtained by intersecting the detector plane and filtering planes. 
Substituting Equation (20) and Equation (23) into Equation (24) and solving for v, we obtain the equation 
of the filtering line on the planar-detector: 

■ufcos s — cos A) 2R 

v(u :s,X) = — - 1 + — . (26) 

sin s sin A sin A 

In the second case, the normal vector to the tangent plane is given by: 

e = (1 — cos s cos A, cos s cos A, cos s sin A). (27) 
The equation of the filtering line on the planar-detector is given by: 

w(cos A — sin s) 2R 

v(u : s, A = '- + — . (28) 

sin A cos s sin A 

4. Weighting Functions of Filtering Lines 

In the following, we define adequate filtering lines for the three-orthogonal scanning case. In order 
to define the different filtering lines, which are necessary to perform the three-orthogonal circular re- 
construction, we first, introduce certain curves, which separate the planar-detector into different regions. 
As an example consider the X-ray source moving on the trajectory C\. We project stereographically the 
trajectories C 2 and C 3 onto the detector plane as shown in Figure 3. Let PC 2 and PC 3 denote these 
projections. The line PC\ is the projection of trajectory C\. Supp f(x) is supposed to be inside a ball 
centered as the origin and with sufficiently small radius r < R, so that the projection of ROI has a circu- 
lar shape, as shown by the shaded region in figure 3. The curves PC 2 and PC 3 split the entire ROI into 
three sub-ROIs: PR m , m = 1,2, 3. The object point x is projected into the area PR m , m = 1,2, 3. Let 
x denotes this projection. 

If i G PRi, L((p) is the projection of the plane through x with normal a x (s,(p). As (p increase, 
a x (s,(p) G P±(s,x) rotates in the clockwise direction on DP(s). From definition of k x (s,(p) with 
equation (14), k is discontinuous if L((p) is parallel to y'(s) or is tangent to C 2 . In the first case, such 
a plane has six points of intersection with C = C\ U C 2 U C 3 , so A; = 1/6 on the side of the jump, 
where y'(s) ■ a > 0 and k = —1/6 on the other side. The filtering line is parallel to y'(s), which is 
denoted as L^see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering 
line Li is w x (s, y?i) = 1/3. In the second case, the number of intersection changes form four to six . 
So the filtering line is tangent to C 2 , which is denoted as L 2 ( see the dot-dashed lines in Figure 3). For 
the signum function in Equation (14) is unchagened, the corresponding weighting function of filtering 
line L 2 is w x (s,(p 2 ) = 1/12. If x E PR 2 , k is discontinuous only L((p) is parallel to y'(s). So the 
filtering line Liis parallel to y'(s). The corresponding weight function is w x (s, (p±) — 1/3. If x e PR3, 
k is discontinuous if L((p) is parallel to y'(s) or is tangent to C 3 . In the first case, such a plane has 
six points of intersection with trajectory, so A; = 1/6 on the side of the jump, where y'(s) • e > 0 and 
k = —1/6 on the other side. The filter line is parallel to y'{s), which is denoted as L x (as shown in 
Figure 3).The corresponding weighting function of filtering line L{\% w x (s,ipi) = 1/3. In the second 
case, the number of intersection changes form four to six. So the filtering line is tangent to C 3 , which is 
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denoted as L 3 ( see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering 
line L 3 is w x (s, cps) = 1/12. Figure 3 summarizes the above information. 



Figure 3. Detector plane with various filtering lines covering projection of ROI shown. 




5. Conclusions 

An exact shift invariant filtered back-projection (FBP) reconstruction algorithm for a cone-beam ver- 
tex trajectory consisting of three- orthogonal to each other circular was derived from Tuy's inversion 
formula. There are several major modifications to Tuy's formula. The first modification is not using the 
Fourier transform with the projection function to deduce the specific inversion formula, but instead using 
the inversion Radon transform. Second, we started with a Tuy-like inversion scheme. Weighting the re- 
dundant data and rewriting the weighted summation into an integral along the source trajectory resulted 
in a shift-invariant FBP reconstruction formula. Last, the "nontangential" condition in Tuy's original 
data sufficiency conditions is relaxed and the nonended condition is further relaxed. In implementing 
of the inversion formula for a cone-beam vertex trajectory consisting of three-orthogonal to each other 
circular, the concrete forms of the filtering lines and according weighting function are the most important 
and difficult segments. In this paper, we obtained above two segments base on analysis about the geo- 
metric properties of projections with the three-orthogonal circular trajectory and the radon planes which 
pass through the reconstruction point. 
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